clc; clear;

% 参数扫描范围
x110_list = linspace(-0.7, 0.2, 40);

% 固定初始条件和参数
q110_0 = 1;
q210_0 = 0;
k = 0.2;
phi0 = 0;
x_init_base = [0; 0; 0; 0; 0];

dt = 0.01;
T_total = 500;
T_trans = 450;  % 暂态时间
steps = round(T_total/dt);
trans_steps = round(T_trans/dt);

% 参数矩阵定义
alpha = 2.4;
A = [-alpha, 0, 0;
      0, 0, 0;
      0, 0, 0];
B = [1, -4, -3.5;
     0, 1, 2;
    -1, -4, 1.5];

n = length(x_init_base) + 4; % 状态空间维度 = 3 + 3 + 3 (这里共9维，但base长度是5，加上4共9维)

LEs_all = zeros(length(x110_list), n);

parfor i = 1:length(x110_list)
    x110 = x110_list(i);
    % 构造完整初始状态向量
    % 顺序：[x110, x120,x130,x210,x220,x230,q110,q210,phi]
    x0 = [x110; x_init_base; q110_0; q210_0; phi0];

    f = @(t,x) mCNN(t, x, A, B, k);

    LEs = computeLEs_QR(f, x0, dt, steps, trans_steps);
    LEs_all(i,:) = sort(LEs,'descend')';  % 转成行向量
end

% 绘制前4个最大LE
figure;
plot(x110_list, LEs_all(:,1), 'r', 'LineWidth',1.5); hold on;
plot(x110_list, LEs_all(:,2), 'b', 'LineWidth',1.5);
plot(x110_list, LEs_all(:,3), 'g', 'LineWidth',1.5);
plot(x110_list, LEs_all(:,4), 'm', 'LineWidth',1.5);
xlabel('x_{110}');
ylabel('李雅普诺夫指数 (LE)');
legend('LE_1','LE_2','LE_3','LE_4');
title('李雅普诺夫指数谱');
grid on;

% ========== 函数部分 ==========

function LEs = computeLEs_QR(f, x0, dt, steps, trans_steps)
    n = length(x0);
    Q = eye(n);
    LEs = zeros(n,1);
    x = x0;

    for i = 1:steps
        % 主轨迹积分 RK4
        k1 = f(0,x);
        k2 = f(0,x + dt/2*k1);
        k3 = f(0,x + dt/2*k2);
        k4 = f(0,x + dt*k3);
        x = x + dt/6*(k1 + 2*k2 + 2*k3 + k4);

        % 计算雅可比矩阵 J(x)
        J = Jacobian(@(z) f(0,z), x, 1e-8);

        % 扰动向量线性演化
        V = J * Q;

        % QR分解，正交化
        [Q,R] = qr(V,0);

        if i > trans_steps
            for j = 1:n
                LEs(j) = LEs(j) + log(abs(R(j,j)));
            end
        end
    end

    LEs = LEs / ((steps - trans_steps)*dt);
end

function J = Jacobian(f, x, eps)
    n = length(x);
    fx = f(x);
    m = length(fx);
    J = zeros(m,n);
    for i = 1:n
        dx = zeros(n,1);
        dx(i) = eps;
        J(:,i) = (f(x + dx) - fx) / eps;
    end
end

function dxdt = mCNN(~, x, A, B, k)
    X1 = x(1:3);
    X2 = x(4:6);
    q11 = x(7);
    q21 = x(8);
    phi = x(9);

    M1 = exp(-q11);
    M2 = exp(-q21);
    Y1 = M1 .* X1;
    Y2 = M2 .* X2;

    W_phi = 1 - 2*exp(phi);
    K = diag([0, 0, k]);

    dX1 = -X1 + A*Y1 + B*X1 + K*W_phi*(X1 - X2);
    dX2 = -X2 + A*Y2 + B*X2 - K*W_phi*(X1 - X2);
    dq11 = X1(1);
    dq21 = X2(1);
    dphi = X1(3) - X2(3);

    dxdt = [dX1; dX2; dq11; dq21; dphi];
end
